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Two-dimensional incompressible viscous driven-cavity flows are computed for Reynolds numbers on the 
range 100-20,000 using a loosely coupled, implicit, second-order centrally-differenced scheme. Mesh se- 
quencing and three-level V'-cycle multigrid error smoothing are incorporated into the symmetric Gauss-Seidel 
time-integration algorithm. Parametrics on the numerical parameters are performed, achieving reductions in 
solution times by more than 60 percent with the full multigrid approach. Details of the circulation patterns 
are investigated in cavities of 2-to-l, 1- to- 1 , and l-to-2 depth to width ratios. 


Nomenclature 

,4 coefficients of linear system 

B f intermediate source terms 

b source term of linear system 

C continuity equation error 

e , e f error vectors 

/ vector of dependent variables 

y source- ter m vector 

g discretized form of g 

h depth of cavity 

/ width of cavity 

P relative pressure 

R f Reynolds number, pVijp 

r, r' residual vectors 

,s arc- length fraction 

s ] grid-clustering parameter 

t time 

U speed of upper (driving) surface 

?’ cartesian velocities 

x , y cartesian coordinates 

ij grid-clustering strength, {3 > 1 

( convergence tolerance 

vorticity 

X s A / / *2 A a* 

A y A t / 2A y 

p viscosity 

p density (constant for incompressible flow) 

a > A / / R e Ax 2 

(T y A / / R e \y 2 

p stream function 

* Aerospace Technologist, Aerothermodynamics Branch, Gas 
Dynamics Division. 


Su bscripts: 

i discrete indicie in ^-direction; i — 1 

j discrete indicie in ^-direction; j = 1, .... J 

c, /, 771 coarse, fine, and medium meshes 
f, x, y differentiation w.r.i. independent variables 

Superscripts: 

k sub-iteration counter 

n time level 

All quantities are non-dimensional. The reference 
values are /, V, pU ", and U/l for lengths, velocities, 
pressure, and time, respectively. When subscripts are 
omitted, they are assumed to be at (?, j). 

Introduction 

V ISCOUS fluid flow in the driven cavity has long 
been a popular test case for evaluating numerical 
techniques. The problem statement is straightforward, 
the geometry simple, yet the governing equations are a 
non-linear system of partial differential equations whose 
stiffness can be adjusted via the Reynolds number, with 
higher Reynolds numbers applying increasing degrees of 
numerical difficulty in achieving a converged steady-state 
solution. 

Recent investigations into this problem by Huser 
and Biringen 1 indicate that by a Reynolds number of 
30,000 the flow has become fully unsteady, u continu- 
ously evolving in time. Shen 2 reports a high-resolution 
scheme that predicts bifurcation of the cavity flow solu- 
tion at Reynolds numbers as low as 10,000. Goodrich ct 
aP also find a periodic solution of this Hopf bifurcation. 
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A partial multigrid scheme has been applied to the 
problem by Nishidaand Nobuyuki" 1 whereby the pressure 
equation was solved with a multigrid procedure and the 
momentum equations are advanced using a Runge-Kutta 
explicit scheme. Another recent solution on the square 
cavity have been put forth by Semeraro 0 using the theta 
scheme. 

The driven-cavity problem statement whose solu- 
tion is sought with the present method was put forward 
by Rubin ci al 6 A loosely-coupled implicit procedure is 
formulated based on symmetric Gauss-Seidel relaxation, 
incorporating three-level multigrid error smoothing on 
each of the governing equations to accelerate conver- 
gence. 

Governing Equations 

Incompressible flow is defined by the continuity 
equation with constant density, 

lly; “t" Vy — 0 ( 1 ) 

in cartesian coordinates. The corresponding perfect-gas 
Navier-Stokes equations' 8 can be written for the prim- 
itive variables, in non-conservative form, as, 

u t + uu T + vuy = -P T + -^-V 2 u (2) 

ti e 

and, 

1 o 

Vt T- uv x -f Wy - -Py + — (3) 
wdiere the Laplacian operator is, 

d 2 d 2 

v '~d^ + w 

Differentiating (Eqns. 2, 3) and incorporating (Eqn. 1) 
yields the Poisson equation for the pressure, following 
Hoffmann/ 1 as, 

V-'/ J = -C, - (u 2 r + 2u y v x + vl) + 2_V 2 C (4) 

where, 

T Uj- T Vy 

is the error of the continuity equation (Eqn. 1), which 
goes to zero in the converged solution. 

Alternatively, the governing equations can be for- 
mulated for the vorticity and stream function. 6 The 
vorticity equation of motion is, 

Ct + U Cr + <y = — ^V 2 C (5) 

and the Poisson equation for the stream function is, 

V“V = C (6) 


The vorticity is related to the velocity components as, 
C = Uy - V x 

and the stream function by, 

Vy = U, VV = — V 


Discretization 

The governing equations (Eqns. 2-6) are dis- 
cretized using a first-order implicit formulation in time 
and second-order expressions for all spatial derivatives. 
Whether applying the governing equations in primitive- 
variable or stream-function/ vorticity formulations, the 
equations take on two basic forms. Setting 

/ = [u, V, C- P, il’] T (7) 

to represent the dependent variables in (Eqns. 2-6), we 
can write for the first three components of (Eqn. 7), 

ft + u fr -f vf y — -tt— V 2 / = —g (8) 

tie 

and for the last two components of (Eqn. 7), 

- V 2 / = -g (9) 

wdiere, 

9 = [Pr, Py , 0, 94 , Cf 

and, 


~g\ - ( \ + u 2 + 2 u y v x + v; - -p-V 2 C 

The discretization strategy will evaluate the com- 
ponents of / implicitly at the advanced timestep, w f hile 
the nonlinear left-hand-side (LHS) terms u and i>, as well 
as the right-hand side (RHS), g , are lagged at the previ- 
ous timestep. creating a loosely-coupled set of equations. 
To reduce notation, / is used to represent f n+1 in the 
difference equations that follow 7 . 

Discretizing (Eqns. 8, 9) yields. 





+ = ( 10 ) 

and, 

_ ^ 1+l — ^ + ^ - Vj + i -2/ + fj-i) 

= -ii (ID 
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The RHS g is the discrete form of y, 


g = 


2A.r 

pn p» 

t i + \ r j- 1 
2 At/ 

0 

</4 

C" 


with, 


C" 

~ 94 ~ _ a7 + 


T+i 




4 Aj“ 
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K+i -«"-i) («"+i - t! f- 1) ( p "+i ~ 

2Aj* Ay 4 Ay 2 
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7 ^ 


rm *ipn , ✓ m rm __ -)fm . f nt 

( i + \ - 4- ( 2 ,j t-j + i + ( J-1 

Aj‘ 2 Ay 2 


and, 


m __ “?>1 


r M = 


C m+1 = 0, 

— I 


i-l ^|_ L j + 1 


Vi 


2A.r ’ 2 Ay 

1’he system of (Eqns. 10, 11) can be rearranged as, 

A ijf + Aj + l fi + 1 + -4 1’ — i / 7 _ i + 


Bi 


( 12 ) 


When solving for », t\ or (,*, the coefficients and RHS of 
(Eqn. 12) are defined as, 

A ij = 1 4- 2(7 T -h 2 cTy 

~ r T A J- U , Aj 1 (Ty it A yl' 


Bij = r - a/ g 


using the notation, 


Ar = 


A / 
2Aj v 

At 


A = 
y 2 Ay 

At 


T ~ R f Ax 2 ' y ~ R e Ay 2 

When solving the Poisson equations for P or V ! i the co- 
efficients and RHS of (Eqn. 12) are defined as, 


4 2 2 
' 4,j -a^ + a7’ 

' 4i±1 = _ A^’ j4j±1=_ A^ 

Finally, the linearized system of governing equations 
(Eqn. 12) is written for interior points as, 


Af = i, 


[2, / — 1] , j = [2, J-l] (13) 


Y 



Upper Surface 
Translates 
at Speed U 


h=l = l 


X 


Fig. 1 Physical domain of driven cavity, 
where the RHS is, 


b = B + fl' 


and B' is used to move the boundary points from the 
LHS to the RHS, 

—Ai-ifi-i, i~ 2 

“^4*+i/t+i < i — 7~ 1 

-Tj-i/y-i, j = 2 

— j + 1 fj + 1 1 / = J — 1 

0, i = [3, /- 2], j = [3, J- 2] 

Thus, (Eqn. 13) is the linear system to be solved at 
each timestep to update the dependent variables from 
time level t = nAt to t = (?i| 1)AE This process is 
repeated until the steady-state solution is achieved. 



Domain 

The domain is a rectangular cavity of width l and 
depth /}, oriented with the lower corner at the origin. 
An infinite upper surface is translated to the right with 
speed U . Figure 1 sketches the physical domain, indicat- 
ing the anticipated flow circulation direction, clockwise, 
for positive V . 

The domain is discretized on a structured cartesian 
grid. Two options for setting the grid-point spacing are 
available; uniform spacing or clustered to walls. The wall- 
clustering function used to generate the arclength frac- 
tion along i=const or j— const lines is, 

i = si (/i+ \)-(i3-l) 

2($i + 1 ) 

where, 


->«-(/ + !) 
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On the left (?’ = I) wail the conditions are expressed 



Fig. 2 Medium-level uniform mesh, 41 by 41 points. 



Fig. 3 Medium-level clustered mesh. 41 bv 41 points. 


as. 


and, 


u " +1 = t - n+1 = o 


P* = 0 => = -pp +1 - -pr + , 

for a second-order extrapolation, or for a first-order ex- 
pression, 

P r = 0 => P n + 1 = P? +] 

For the stream function and vorticity. 


iA»+ 1 


0 


vv = o 


— V T 


Cl = C>i 


C 


n+1 


2 


Ax 


2 V t'+l 


Similarly, on the right (?*=/) wall, 
u rl + 1 = i ,n+1 = 0 

Pr = 0 => / ,T ‘ + 1 = Jp/Ll ~ 1/^2 

O <J 


or, P n+] = pr\ 


+1 


! fi + 1 


V" ■ - = 0 

0 = 0 => t/’I'+i = 1 

2 


C = t'xx => C n+1 = 


On the bottom ( j = 1) wall, 


Ax 2 


Vi - 1 


The typical three-level multigrid procedure employed 
here solves the governing equations on an 81 by 81 fine- 
mesh and smoothes the error on two additional meshes, 
dimensioned 41 by 41 and 21 by 21, generated from 
the fine mesh by recursively removing every-other point. 
Sample 41 by 41 medium-coarseness meshes, which re- 
produce clearer than the fine meshes, are presented in 
Figure 2 for a uniform grid and Figure 3 for a grid clus- 
tered with (3 = 1.6. The 81 by 81 fine mesh can be 
inferred from Figures 2 and 3 by simply doubling the 
number of grid lines in both i and j directions. 

Boundary and Initial Conditions 

Viscous, solid walls are imposed on all four sides of 
the cavity. This imposes a no-slip condition on the ve- 
locity components at the wall. For pressure, the wall- 
normal pressure gradient is set to zero. This is the 
boundary-layer assumption, which becomes more accu- 
rate as the Reynolds number increases. The stream func- 
tion is set to zero on the boundary. 


ti n+1 = t;” +1 = 0 


Py = 0 


1 - 


j+2 


or, P' , + 1 = p;* +1 


0 n+1 = 0 


= 0 

C z Vyy 


,!, n — ...” 

~ »J + 1 






jV j+ i 


l he top ( j = J) wall is translating horizontally, so, 

u n+l = V, C ‘ +1 = 0 


Py= 0 


=> P n+l = - \pf - 2 


or 


v y = V 


pn + l _ pn_ i 
V ;,1 + 1 = 0 

* V ’"+1 = 2 1 'Ay + 
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C = 4'yy => C + ' = ^(U^y+fj-l) 

The boundary conditions are updated explicitly at 
each timestep. 

The flowfield is initially at rest at 1 = 0”, with all 
dependent variables set to zero, except for the pressure 
which is set to one everywhere. The upper boundary is 
then impulsively started at t = 0+ to its full speed, U . 
Thus, the solution is started with the upper surface set 
to u — U or C = 2 U/Ay. 

Solution Procedure 

Evolution in Time 

Beginning with the initial conditions at time i ~ 0 + , 
the penta-diagonal linear system (Eqn. 13) is approxi- 
mately inverted at each time step to advance the depen- 
dent variables to their steady-state values. This inver- 
sion is performed in a loosely-coupled fashion by iterat- 
ing for the first, second, and then fourth components of 
(Eqn. 7) at each timestep when solving the primitive- 
variables formulation, or iterating for the third and then 
fifth components of (Eqn. 7) when solving the stream- 
function/ vorticity formulation. Coupling between the 
equations is achieved when the A matrix and b vector of 
(Eqn. 13) are reformulated at each new time level. 

Local time stepping is employed by setting At at 
each point to its maximum allowable value based on the 
input constraints on the two stability parameters A maJ 
and cr m ar ? such that, 

x _ . ( 2 A ma j. - m i n ( A x ? y , A yij ) \ 

- nun ^ . min( (Aatjj) 2 , (At/ij) 2 ) ) 

Local time stepping reduces the time accuracy of the 
solution on non-uniform meshes but is commonly em- 
ployed to speed convergence to a steady-state result. 
The present algorithm already suffers poor time accu- 
racy due to the first-order temporal derivative, lagged 
coefficients, loose coupling between governing equations, 
and approximate inversion of the linear system at each 
timestep. However, it is the steady-state solution that 
the present method seeks, whose accuracy is determined 
by the second-order formulation of the spatial deriva- 
tives, and not the time-history of the flow. 

Relaxation Algorithm 

A symmetric Gauss-Seidel (SGS) iteration sweep is 
performed to approximately invert Af — b. Forward and 
backward sub-iteration sweeps are performed for each 
component of / before moving on to the next component. 
This involves f \ , f->, and for the primitive- variable for- 
mulation, or /a and / 5 for the stream-function/vorticity 
formulation. The forward sweep loops on i for each j, 


while the backward sweep loops on j for each /, attempt- 
ing to minimize any bias in the solution from the SGS 
sweep directions. 

A forward half-sweep of the SGS procedure on 
(Eqn. 12) at the k th sub-iteration is, 



The backward portion of the symmetric sweep follows 
the same pattern as the forward sweep, but using the k th 
values of the i + 1 and j + 1 variables. Typically, good 
performance is achieved with only 1-4 sub-iterations. 

Grid Sequencing 

Grid sequencing, or nested iterations, is employed 
at the beginning of the solution to speed convergence 
on the finest grid. The initial, coarsest grid is obtained 
from the fine mesh by retaining only every fourth point 
in both the i and j directions, for a factor of 16 reduction 
in the total number of grid points. The initial conditions 
are advanced on this grid for a number of timesteps equal 
to I c J c / 2. 

The coarse- mesh solution is then prolongated to the 
medium mesh, which retains every-other grid point from 
the fine mesh. The prolongation employed here is a di- 
rect injection for overlapping points and two-point or 
four-point averaging for the interstitial points. Four 
times as many timesteps are then taken on the medium 
mesh as on the coarse mesh, and the result ing solution is 
then prolongated to the fine mesh, completing the grid 
sequencing. 

Multigrid Cycle 

A full three-level multigrid V-cycle is employed to 
damp low frequency errors in the solution. SGS relax- 
ation on the fine mesh effectively damps the high fre- 
quency error content, but poorly damps low frequency 
errors. 10 The multigrid cycle damps these low frequency 
errors by restricting the residual to progressively coarser 
meshes, aliasing the low frequencies to higher frequen- 
cies which are then damped effectively by the SGS re- 
laxation algorithm. The error is then prolongated back 
to the fine- mesh solution as a correction to the depen- 
dent variables. 

The three-level V-cycle begins by evolving the solu- 
tion in time on the finest grid for a number of t imesteps, 
typically 5-10. Then the residual is formed as. 

rj=b-Af 

The residual is restricted, via simple point-to-point ex- 
traction, to the medium mesh, where the error equation 
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is relaxed, 


Ae m — r m (14) 

A new residual is formed for this error equation as, 

~ r rji — Ae nt 

and is restricted to the coarse mesh, where another error 
equation is relaxed, 

Ae c = r' 

Once the error on the coarse mesh has been 
smoothed through several SGS sweeps the error e c is 
prolongated to e f m on the medium mesh, and is applied 
as a correction to e m , 

* £' m -b e m 

Again the error is smoothed on the medium mesh 
(Eqn, 14) for several SGS sweeps and then the error e m 
is prolongated back to e'j on the fine mesh. The solution 
is now corrected by this error, 

f-f + e'j 

and is evolved in time again on the fine mesh for several 
timesteps. The multigrid cycle is then repeated until 
convergence is achieved for the fine-mesh solution. 

Convergence Criteria 

The ^ 2 -norm of the change in the dependent vector, 
/, between time-levels n and rc + 1 is monitored to judge 
convergence, 

nr +1 -/ n ii 2 <* 

Lo-norms are normalized by their value at the first 
timestep, with a typical value for c of 10“ 4 . 

Results 

Test Case 

A standard test case is chosen to compare the com- 
putational performance of the various techniques em- 
ployed in the present method. Table 1 lists the param- 
eters which are desired to be optimized using the test 
case. 

The test case is chosen on the unit cavity, h — l = 1, 
with a uniform 81 by 81 fine mesh. The upper surface is 
translated to the right at speed U — i, and the Reynolds 
number is 1000. The stream-function/vorticity form of 
the governing equations are solved. 

Straight SGS relaxation without mesh sequencing 
or multigrid was used to generate the baseline solutions 
and optimize the maximum allowable timestep. For this 
case the timestep limit on At from A maT is about six 


Table 1 Algorithmic parameters to optimize. 


Parameter 

Test Values 

A mari &max 

1, 2, 3, 4 

SGS sub- iterations 

1, 2 

Mesh sequencing 

none, max(/, ./), I+J, 
/ -.7/2, / 7/4, / -7/8 

V -cycle iteration pattern 

none, 3/-2 m -l f -2 m -3/, 
6f-4 m -3 c -4 m -Q f: 
2/-3 m -4 r -3 m -2/ 

Grid refinement 

41x41, 81 x 81 



Fig. 4 ( 'divergence rates for A maj ., <T mar variations. 

times more restrictive than from <r mar . Figure 4 com- 
pares convergence rates for various values of X max and 
o’majr. With A m ax — 4 the scheme is unstable, diverging 
after 50 timesteps. With a A mar = 3, the solution does 
not diverge, but is clearly on the ragged edge of stabil- 
ity. Converging three orders of magnitude in 357 CPU 
seconds, A T1iaT — 2 gives the best convergence rate while 
maintaining stability, and is used for the remainder of 
the study. For smaller values of A max the convergence 
rate slows. Convergence rates, measured in CPU sec- 
onds, were obtained on a 40 MHz SUN Sparcstation 2 
for the present study. Single-precision numerics, approx- 
imately seven significant digits, was used. 

Figure 5 compares the effect of performing multiple 
SGS sweeps during each timestep. Performing two sets 
of forward and backward sweeps during each SGS cycle 
is seen to take longer to reach the converged steady- 
state solution, 400 versus 357 CPU seconds, than tak- 
ing single forward-backward sweeps in each SGS cycle. 
With two SGS sweeps, A majr = 2 proved to be unstable, 
so A max = 1.5 was used. This correspondingly smaller 
timestep is contributing to the two-pass SGS relaxation 
taking longer to converge. 
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Fig. 5 ( "divergence using one and two SGS 
sub-iterations. 



Fig. 6 Convergence rates with mesh sequencing. 

Adding mesh sequencing to the solution procedure 
reduced the convergence time in half, to 147 CPU sec- 
onds. Several methods for specifying the number of 
coarse-grid iterations to perform were tested, and are 
listed in Table 1. The corresponding convergence rates 
are presented in Figure 6, where the general trend is that 
the more iterations taken on the coarser grids, the faster 
the overall convergence. Setting the convergence limit as 
i — 10 -3 , the fastest convergence is obtained using, 

Ir 'Jr _ 1 -J 

4 “ 64 

steps on the coarse grid. The La-norm jumps when the 
solution is prolongated to the medium mesh, where, 

Irn * Jrn = / J 

4 16 

steps were taken before prolongating to the fine mesh. 
Similar trends are seen for the other functions used to 



Fig. 7 Mesh sequencing convergence rates. 

determine the number of steps to take on the nested 
grids. 

While Figure 6 indicates the 1 J/4 function gives 
the fastest convergence in CPU time, the function / -J / 2 
has a steeper convergence slope at e — 10~ 3 . Carrying 
these two functional variations to further convergence. 
Figure 7, shows that the functional form 7 -J/2 converges 
the fastest toward t — 10” 4 , and is the chosen function 
for the present method. 

While mesh sequencing was found to speed con- 
vergence on the fine mesh, Figure 7 shows that the 
asymptotic convergence rate to small ( is st ill fairly slow. 
Adding multigrid to the solution procedure increases the 
rate of convergence to small ( by more effectively damp- 
ing all frequency contents in the error. This change in the 
slope of convergence toward small ( using the multigrid 
approach is demonstrated in Figure 8, where all three 
variations in l -cycles for the multigrid implementation 
converge to c — 10~ 4 in 450-500 CPU seconds, while 
the non-multigrid solution has reached an asymptotic 
convergence rate that will take at least twice as long to 
converge to ( — 10 -4 . The convergence histories of the 
multigrid solutions display spikes each time the error is 
prolongated back to the fine grid as a correction. While 
these correction steps temporarily increase the L 2 -norm, 
they clearly increase the global rate of convergence. All 
three V -cycles tested performed nearly the same, with 
the 3/-2 m -l r -2 m -3f pattern performing the best, con- 
verging four orders of magnitude in 466 CPU seconds, 
and is the pattern chosen for the present method. 

The velocity vectors for this cavity flow problem, 
as computed with the fastest-converging multigrid so- 
lution, are drawn in Figure 9. For clarity, only every 
fourth vector in both i and j directions has been plot- 
ted. The highest velocities are seen on the upper, driving 
surface, as expected. A strong circulation region forms 
on the interior of the cavity at this Reynolds number, 
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Fig. 8 Multigrid convergence rates with variations in 
V -cycle. 


Y 



R e = 1 0 0 C 

81x81 uniform 
mesh 

Plot only every 
fourth point . 


X 


Fig. 10 Streamlines in viscous driven cavity. 
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Fig. 11 Velocity component profiles through the core. 


Fig. 9 Velocity vectors in viscous driven cavity. 

R ( = 10 3 . The circulation resembles are solid-core ro- 
tation, with the highest velocities on the perimeter and 
decreasing velocities toward the center, as opposed to 
a vortex, which would have increasing velocities toward 
the eye. The bottom corners are regions of relatively 
stagnant flow, and the upper left corner shows the suc- 
tion characteristic of the entrainment into the upper, 
driven boundary layer. 

Figure 10 traces the corresponding streamlines. The 
primary circulation is seen to encompass nearly the en- 
tire cavity. Two recirculation regions are formed in each 
of the lower corners, but, recalling Figure 9, the velocity 
magnitudes are small in these recirculation regions. 

Velocity components through the boundary layers 
are extracted in figure 11. The u component has been 
extracted along a vert ical slice through the center of the 
domain while the v component is taken along a hor- 
izontal slice, again through the center of the domain. 
The central circulation region is confirmed to resemble a 


solid-core rotation, with linear velocity profiles through 
the core. Considering the grid-point distribution, the in- 
terior of the domain is well resolved. The bottom and 
left-wall boundary layers each have 15 points through the 
layer, providing moderate resolution. The upper surface 
and right wall have only 8-9 points through the laminar 
boundary layer. 

Grid Convergence 

The test case is repeated, on a uniform 41 by 41 
mesh. Velocity components through the cavity are plot- 
ted in Figure 12, along with the 81 by 81-mesh solution 
from Figure 11. T he 41 by 41 mesh is clearly not re- 
solving the boundary-layer gradients nor the maximum 
velocities on the periphery of the circulation region. The 
linear trend in velocity variation through the circulation 
is picked up on the coarser grid. Since the 41 by 41-mesh 
solution does not generate the correct maximum veloci- 
ties in the circulation, nor have the correct slope of t he 
velocity profiles near the walls, and hence an inaccurate 
skin friction, the solution is not grid converged for this 
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Fig. 12 Coarse- grid solution to viscous- cavity test 
case. 



Fig. 13 Grid convergence for driven-cavity flow, 
mesh. 

The test case is repeated again, now on a 121 by 121 
mesh. Figure 13 over-plots this very-fine mesh on the 
81 by 81 solution. Excellent agreement is seen in the 
velocity profiles across the cavity. The boundary-layer 
profiles show good agreement as well on all walls, sug- 
gesting the 81 by 81-point mesh sufficiently resolves the 
flow at this Reynolds number. 

Grid Clustering 

The presence of the high velocity gradients in the 
boundary layers, as compared to the gradients on the 
interior of the cavity, suggests the use of a stretched grid 
to resolve the flow with fewer points, and hence in less 
time and computer resources. 

Figure 14 over-plots the u and v velocity profiles 
across the cavity for the 81 by 81 uniform grid and 
a 61 by 61 clustered mesh. Both solutions are at a 
Reynolds number of 1000, and the clustered grid was 



Fig. 14 Clustered-grid velocity profiles. 

generated with j3 — 1.6. The basic flow pattern remains 
the same with the clustered grid, with good agreement 
through the upper boundary layer and across the pri- 
mary circulation zone. However, the agreement is not as 
good in the bottom boundary layer and at the periphery 
of the circulation region, where the velocity peaks are 
under- predicted. 

Better agreement is sought in the grid-clustered so- 
lutions. but the general trend of a degradation in the so- 
lution in the presence of stretching between mesh lines 
was seen in other cases as well. Sources of this error 
could arise from the single- precision numerics, the for- 
mulation of the algorithm for solving t he governing equa- 
tions, or in an unlocated programming bug. The inabil- 
ity to achieve reliable solutions with grid-point, cluster- 
ing to the walls limits the practical Reynolds number to 
which the present method can be applied. 

High- Reynolds- Number Flows 

Driven flow in the unit cavity at, a Reynolds num- 
ber of 2000 was computed on a 101 by 101 mesh. Cavity 
streamlines are drawn for this case in Figure 15. The 
streamline pattern for this case remains materially the 
same as for the R e — 1000 solution in Figure 10. Velocity 
profiles for R e = 1000 and R e = 2000 are plotted in Fig- 
ure 16 at every fourth (/, j) point. The linearity in the 
velocity profiles through the circulation region is begin- 
ning to be lost at the higher Reynolds number and the 
boundary layers have thinned, but overall the flow has 
remained much the same with a doubling in Reynolds 
number. 

Further increasing the Reynolds number to 5000 
produces a dramatic shift in the flowfield. Figure 17 
plots Velocity vectors for this case. Clearly, the pri- 
mary central circulation has broken down, and is re- 
placed by a clockwise circulation in the upper right cor- 
ner and a counterclockwise recirculation of comparable 
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Fig. 16 Velocity profiles at higher Reynolds number. 


Fig. 19 Velocity vectors at R< — 10,000. 
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Fig. 17 Velocity vectors at R P = 5000. 

size, though of lower velocity magnitudes, near the upper 
left corner. The corresponding streamlines are traced in 
Figure 18, where the dual-circulation pattern is again 
observed. 


Solutions were obtained as well for R e = 10, 000 and 
R e = 20,000, the highest Reynolds number for which a 
reliable solution was obtainable in a reasonable amount 
of time, about one hour, with the present method. Veloc- 
ity vectors for these two cases are plotted in Figures 19 
and 20 at every-other (?*, j) point. The circulation pat- 
tern seen at lower Reynolds numbers is breaking down, 
and the magnitude of velocities on the interior of the 
cavity are decreasing. The boundary layers are thinning 
as well, and are poorly resolved on the 101 by 101 mesh 
at R e = 20,000. 


The corresponding streamlines are traced in Fig- 
ures 21 and 22. The primary circulation zone lias contin- 
ued to shrink in size in going from Reynolds numbers of 
1000, to 5000, 10,000, and 20,000. Also, the counter- 
rotating recirculation has moved, from the upper-left 
corner at R f = 5000, to a more central location at 
R e = 10,000, and finally all the way to the upper-right 
corner by R f — 20.000. 
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Fig. 20 Velocity vectors at R e — 20,000. 


Fig. 23 Tail-cavity flows. 
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Fig. 24 Shallow-cavity velocity vectors. 


Fig. 21 Streamlines at R. e — 10.000. 
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Fig. 22 Streamlines at R f — 20,000. 


Effect of Cavity Aspect Ratio 

The cavity dimensions were changed so that the 
depth is twice the width. Solutions were generated at 
Reynolds numbers of 100 and 1000 on uniform 41 by 61 
grids. Velocity vectors for these two cases are plotted 
in Figure 23. The circulation core is closer to the driv- 
ing surface in the more viscous, ^=100, solution, while 
the higher Reynolds-number solution exhibits generally 
higher velocities around the circulation region. 

The flowfield changes character when the cavity is 
made twice as wide as deep. Figure 24 plots velocity vec- 
tors, at every fourth (/, j) point, in a cavity with h = 1 
and / = 2 at a Reynolds number of 1000. An 81 by 61 
uniform mesh was used for this calculation. The strong 
primary circulat ion forms with a nearly circular shape at 
the right side of the cavity. The recirculation region on 
the left side of the cavity is approximately half the size 
of the primary circulation, though with much reduced 
velocities. Streamlines for this case are traced in Fig- 
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Fig. 25 Shallow-cavitv s t ream 1 i n es . 

ure 25, where a second, weak recirculation can be seen 
in the lower right corner. 


Summary of Results 

Incompressible viscous-flow solutions have been 
computed for the two-dimensional driven-cavity prob- 
lem. The governing equations for perfect-fluid viscous 
flows have been discretized in a finite-difference for- 
mulation in a loosely coupled implicit scheme, formu- 
lated either in terms of the primitive variables or the 
stream-function and vorticity.. Mesh sequencing and 
full, three-level, multigrid error smoothing are incor- 
porated to speed convergence of the symmetric Gauss- 
Seidel time-integration algorithm. 

Convergence was found to be faster when using 
a single, rather than multiple, symmetric Gauss-Seidel 
sweep at each time level to approximately invert the lin- 
ear system. Several mesh sequencing iteration functions 
were tried, and the best was able to reduce solution time 
by 58 percent over an un-sequenced solution. The addi- 
tion of the multigrid procedure further accelerated con- 
vergence on fine meshes. 

At a Reynolds number of 1000 the cavity flow is 
characterized by a large solid-core circulation, encom- 
passing nearly the entire domain. Counter-rotating cir- 
culations are found in the lower corners of the cavity, 
but with small velocities. As the Reynolds number is 
increased, the primary circulation is driven toward the 
upper-right corner, with the recirculations moving up 
and across the cavity. The highest Reynolds-number 
flow computed was at 20,000. 

Vertical stretching of the cavity was not found to 
alter the flow patterns significantly, while horizontal 
stretching of the cavity created a large recirculation re- 
gion, at a Reynolds number of 1000, and a laterally 
stretched entrainment pattern. 
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